Finite differencing second order systems 
describing black hole spacetimes 
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Keeping Einstein's equations in second order form can be appealing for computational efficiency, 
because of the reduced number of variables and constraints. Stability issues emerge, however, which 
are not present in first order formulations. We show that a standard discretization of the second 
order "shifted" wave equation leads to an unstable semi-discrete scheme if the shift parameter is 
too large. This implies that discretizations obtained using integrators such as Runge-Kutta, Crank- 
Nicholson, leap-frog are unstable for any fixed value of the Courant factor. We argue that this 
situation arises in numerical relativity, particularly in simulations of spacetimes containing black 
holes, and discuss several ways of circumventing this problem. We find that the first order reduction 
in time based on "ADM" type variables is very effective. 

PACS numbers: 02.70.Bf, 04.25.Dm 



I. INTRODUCTION 

In recent years there has been a growing interest in 
discretizing the second order Einstein's equations, in the 
harmonic gauge or its generalization,without reducing 
the system to first order form Q, 0, Q ■ The reduction 
process requires the introduction of auxiliary variables 
approximating first derivatives of the fields and the in- 
troduction of additional constraints. Whereas there are 
clear advantages in keeping the system of equations in 
second order form, including the fact that local well- 
posedness of the continuum Cauchy problem has been 
shown and the expectation that in general this would 
lead to smaller numerical errors we point out diffi- 
culties that can arise when a standard discretization is 
used. 

After analyzing a toy model problem that captures the 
essential difficulty, and pointing out its relevance to nu- 
merical relativity, we discuss different solutions to this 
problem. The first order reduction in time based on the 
introduction of "ADM" type variables seems to be the 
most attractive of these solutions. 



II. THE SHIFTED WAVE EQUATION 

We start with the wave equation in one spatial dimen- 
sion, <f>n = 4>xx, and perform a Galilean change of coor- 
dinates, t = t , x = x — (3t, where /3 is a constant. This 
leads to 



kt = 20<h* + (1 - /? 2 )<^ 



(1) 



which wc will refer to as the shifted wave equation. By 
performing a differential reduction to first order we see 
that the characteristic variables and speeds are </>(— /3^> x ± 
4> X ) /3i 1. The variable (f> is also a characteristic variable, 
the speed of which is undetermined (it depends on the 
details of the reduction and one can choose what one 
pleases). The initial value problem for this system is 



well-posed for any value of (3. In fact, an energy estimate 
can be obtained by noting that the quantity 



b x f + <&) dx 



(2) 



is positive definite in <fi t , <f> x and is conserved for any (3. 

We introduce the grid Xj = jh, where h is the space 
step, and the grid-function 4>j{t) approximating <fi(t,Xj). 
Leaving time continuous, the standard second order ac- 
curate approximation of Eq. is 



dt 2 



(3) 



where hD + Uj = 



u 7 -_i and 



2Dq = D + + D- . Consider the discrete quantity 

{<fH,<h)h + O--0 > KD+4>,D+<l,) h , (4) 

where (it, v)h = J2j u j v j^ 1 - For |/3| < 1 this expression 
is positive definite in (j> t , D + <j> and is conserved j^. As 
in the continuum case, the energy estimate follows. The 
semi-discrete system is stable. 

On the other hand, if \(3\ > 1, there does not exist 
a positive definite quantity from which one can derive 
a discrete energy estimate. A closer look at Eq. (0 re- 
veals that there might be a problem with the highest 
frequency (and those nearby), due to the fact that Dq 
is unable to see it, Do(—iy = 0. Consequently, at this 
frequency Eq. @ appears to be elliptic. It is not difficult 
to show that the semi-discrete problem admits solutions 
that grow exponentially without bound in h. Inserting 
4>j{t) = e st cf>j into Eq. © we obtain 



s 2 



<t>3 = SPi&j+l ~ <Aj-l) + (1 - P^i^J + l - 2( t>j + 4>j-l) 



where s = sh. For <j>j = (— l)- 7 , we get s 2 = A{(3 2 — 1). 
Hence, the grid-function 



(5) 



2 



is a solution of Eq. @, the growth of which cannot be 
bounded independently of h. Notice that this analysis 
also applies to the first order in time, second order in 
space system 



dt 
dt 



T. 



2pD T 3 + {l~fi 2 )D + D_ 



(6) 
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In particular, this shows that schemes such as the forward 
Euler, backward Euler, Runge-Kutta, Crank-Nicholson 
and leap-frog methods applied to either Eq. @ or system 
(JSJ) arc unstable if \(3\ > 1. The scheme is also unstable for 
|/3| = 1. However, in this case the instability is less severe 
(the system admits linearly growing frequency dependent 
solutions). 

A toy model problem similar to the shifted wave equa- 
tion was considered by Alcubicrre and Schutz 0, who 
proved instability for an implicit scheme and proposed 
using causal differencing @, H, El EH to eliminate the 
instability. Our semi-discrete analysis leads to a more 
general result, namely that the instability is due to the 
spatial discretization and does not depend on the time 
integration. Furthermore, it is important to realize that 
this type of instability does not appear in fully first or- 
der systems. In these cases one can handle high values of 
the characteristic speeds by choosing a sufficiently small 
Courant factor (and possibly adding artificial dissipation 
in the variable coefficient case). Whenever causal differ- 
encing has been applied to first order systems, it has not 
brought substantial improvements 

Before discussing how we propose to fix this problem, 
we show how it can arise in discretizations of second or- 
der systems of Einstein's equations. For concrcteness, we 
consider formulations having principal part determined 
by a wave operator of the form g^d^dv, where g^ v is 
the inverse 4-metric of spacetime. Precisely this oper- 
ator appears in the (generalized) harmonic gauge 
We keep the system in second order form and use the 
standard spatial discretization. Assuming that the co- 
ordinates are chosen such that the t = const, slices are 
space-like, i.e., g u < 0, one can expect the instability to 
arise whenever the spatial coordinates are such that an 
x l = const, hyper-surface is also space-like, i.e., g il < 
(no sum). Again, to the highest grid frequency this prob- 
lem appears to be elliptic. Interestingly, the last condi- 
tion, g 11 < 0, is a requirement for excision, a technique 
often used in numerical relativity to eliminate the black 
hole singularity from the computational domain. This 
shows that when discretizing second order systems de- 
scribing spacetimcs containing black holes, one has to 
ponder over the discretization. 

Another instance in which this type of instability can 
arise is when rigidly co-rotating coordinates are used. 
These coordinates are introduced to attempt to keep a 
binary black hole system at a fixed coordinate location 
[l2L . At large distances the semi-discrete wave opera- 
tor effectively becomes elliptic for the highest frequencies. 



We now go back to the model equation (J3J) and discuss 
several methods to overcome the instability that occurs 
for |/3| > 1, without reducing the spatial derivative. 

Method 1: We know that the addition of artificial 
dissipation can sometimes stabilize otherwise unstable 
schemes. If we modify system (J5J) as follows 



(7) 



ah 3 (D+D-) 2 T, 



j > 



= Tj - o-h 3 (D + D-) 2 (j)j , 

= 2(3D Tj + (1 - [3 2 )D + D^ 

we see that the von Neumann condition, which is only 
a necessary condition for stability, is satisfied for suf- 
ficiently large values of the dissipation parameter (for 
example a > 0.385 for |/3| = 2). However, such a 
great amount of dissipation demands high resolution to 
prevent excessive damping and requires a rather small 
Courant factor. For fourth order Runge-Kutta (4RK) in 
the |/3| = 2 case we need k/h < 0.289, where k is the 
time step. 

Method 2: Perhaps the simplest solution is to replace 
the one sided operators D± in Eq. with the centered 
one, Dq. This amounts to discretizing the second spatial 
derivatives with the D 2 , operator instead of D+D-, as 
suggested in 0, ^J, leading to a scheme with a five 
point stencil instead of three. With such discretization 
the discrete version of @, with the replacement d x — > 
Do, is conserved and a von Neumann stability analysis 
gives a Courant limit of ^8/(1 + |/3|) for 4RK. 

At first glance this method appears to be very effective. 
It suppresses the exponentially growing mode (JSJ and it 
allows for a rather large time step. However, the fact 
that Do is blind to the highest frequency means that the 
discrete conserved quantity is unable to capture it and, 
as discussed in greater detail in 0], the method is not 
robust in the sense that a perturbation of the equation 
by lower order terms can trigger (exponentially growing) 
numerical instabilities. Although it is possible that arti- 
ficial dissipation may cure this problem, this needs to be 
explored. Whatever the case may be, a five point stencil 
is likely to unduly complicate the treatment of bound- 
aries. 

Method 3: Another alternative is to rewrite Eq. JQ) as 
d t <j> = pd x <j> + Ii, (8) 

d t n = f3d x n + d*ct), 

where we have introduced the variable n = dt<t> — f3d x (j). 
The standard second order accurate discretization now 
gives 



dt-^ 



n 



3 1 



(9) 



--/?A)Jn, = D+D-fa. 

Note that in terms of the original second order system 
this spatial discretization corresponds to 

^ - 2/3A)^ + P 2 D 2 ^ = D + D-4>j , (10) 
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which has a five point stencil. Incidentally, for large (3 
it is not possible to construct a centered, second order 
accurate, three point stencil, stable approximation of the 
second order equation Q , without performing a first or- 
der reduction in time. System is stable for any value 
of /?, as it conserves the discrete quantity 

{n,u) h + {D + <j>,D + 4>) h . (ii) 

With 4RK and for large (3 it has a Courant limit compa- 
rable to that of method 2. In particular, for \(3\ = 2 we 
get the condition k/h < 0.803. Furthermore, the numer- 
ical speeds of propagation associated with system (JHJ arc 
closer to the exact ones than those of method 2. 

Both the continuum system and approximation (j5J arc 
non-dissipative. They admit a conserved energy. If the 
advective terms (the terms multiplied by (3) in the semi- 
discrete system (|§J) are approximated with second or- 
der accurate one-sided operators D + (l — |-D+) rather 
than D , assuming f3 > 1, wc are trading a conservative 
scheme with one which is dissipativc. Although in the 
variable coefficient case this may be effective in obtain- 
ing stability, in this particular case, with 4RK and (3 = 2, 
the scheme requires k/h < 0.332, which is less than half 
what is needed by the centered approximation. 

Finally, we point out that in the fourth order accu- 
rate case, d x -> D (l - \h 2 D + DJ), d\ -> D+D-(l - 
j^h 2 D + D_), the results of this paper remain qualita- 
tively unchanged. For > 1 the standard discretiza- 
tion of Q is unstable, unless copious amount of artifi- 
cial dissipation is added to the scheme[2lj, whereas the 
discretization of © is stable for any value of the shift 
parameter. 

III. CONCLUSION 

Our analysis demonstrates that when using formula- 
tions which are second order in space one has to exercise 



caution, even in the absence of boundaries. We find that 
the instability discussed in Q , which motivated the intro- 
duction of causal differencing, is due to the mixing of the 
Do and D + D_ operators in the spatial discretization and 
therefore only appears with second order in space sys- 
tems. The fact that it arises already at the semi-discrete 
level, as in Eq. ©, shows that no time integrator can 
fix it, not even implicit ones. One can expect such insta- 
bility to arise in numerical relativity simulations based 
on standard spatial discretization of fully second order 
systems near black holes and at large distances from the 
center of a rigidly rotating coordinate system. 

We believe that a simple and effective method of elim- 
inating the instability consists in rewriting the system in 
"ADM" form before discretizing it, as in ©. When this 
is done, the resulting approximation, which is still cen- 
tered, is stable for any value of the parameter (3. It is 
straightforward to prove stability using the discrete en- 
ergy method and higher order accurate schemes can be 
easily constructed. Interestingly, the structure of system 
ijHJ is similar to commonly used first order in time, sec- 
ond order in space formulations of Einstein's equations, 
such as the Baumgartc-Shapiro-Shibata-Nakamura sys- 
tem [num. 
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